In [1]:
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib as mpl
import matplotlib.cm as cm
mpl.rcParams['figure.figsize'] = (7, 7)
plt.style.use('bmh')
mpl.__version__
import skimage
from skimage import filter
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage.draw import circle_perimeter
import glob
import numpy as np
In [2]:
from saaspy.image import image, find_center
from saaspy.collection import collection
In [3]:
dir = "/Users/schriste/Desktop/SAAS/"
list of files taken during calibration. See SAAS labnotebook for a list of what each one is.
In [4]:
fits_files1 = glob.glob('/Users/schriste/Desktop/SAAS/FOXSI_SAAS_141109_*.fits')
fits_files2 = glob.glob('/Users/schriste/Desktop/SAAS/FOXSI_SAAS_141110_*.fits')
fits_files = fits_files1 + fits_files2
In [5]:
fits_files
Out[5]:
In [6]:
cl = collection(fits_files[2:])
In [7]:
cl.data
Out[7]:
Checking out one random file
In [8]:
index = -10
s1 = image(fits_files[index])
print(fits_files[index])
In [9]:
s1.imshow()
s1.overlay()
plt.colorbar()
Out[9]:
In [10]:
plt.imshow(s1.data, origin='upper', cmap=cm.Greys_r)
plt.colorbar()
Out[10]:
In [11]:
s1.hist()
In [12]:
s1.roi_auto()
s1.imshow()
plt.colorbar()
Out[12]:
In [13]:
s1.roi
Out[13]:
In [14]:
#find_center(s1, sigma=0.75)
Gathering all of image maxima
In [15]:
cl.data.plot(subplots=True)
Out[15]:
From the plots below you can see when the position was stable (from index 5 to 13) and when the offpointing occured. This is corraborated by comparing to the X-ray alignment images
In [16]:
x = cl.data['mindex_x'].values
y = cl.data['mindex_y'].values
fig, (ax1, ax2) = plt.subplots(ncols=2)
ax1.plot(x)
ax2.plot(y)
plt.show()
using only image max
In [17]:
print(np.average(x[5:13]), np.std(x[5:13]))
In [18]:
print(np.average(y[5:13]), np.std(y[5:13]))
The camera was moved during staking which invalidates this value
definitely need some tweeking
In [19]:
#x_hough = []
#y_hough = []
#i = 0
#for f in fits_files[2:]:
# s = saas_img(f)
## x_hough.append(s.find_center(sigma=0.8)[0])
# y_hough.append(s.find_center(sigma=0.8)[1])
# print(i, x_hough, y_hough, f)
# i = i + 1
In [20]:
s = []
plate_scale_files = ['FOXSI_SAAS_141110_003437.fits',
'FOXSI_SAAS_141110_003449.fits',
'FOXSI_SAAS_141110_003516.fits']
for f in plate_scale_files:
s.append(image(dir + f))
In [21]:
s[0].imshow()
In [22]:
s[1].imshow()
In [23]:
s[2].imshow()
Probably the best image
In [24]:
best_img = s[0]
best_img.imshow()
In [25]:
best_img.roi = [500, 1000, 400, 800]
best_img.imshow()
x1 = [752, 914]
y1 = [756, 438]
x2 = [752, 928]
y2 = [756, 423]
plt.plot(x1, y1)
plt.plot(x2, y2)
Out[25]:
In [26]:
r1 = np.sqrt((x1[0] - x1[1]) ** 2 + (y1[0] - y1[1]) ** 2)
r2 = np.sqrt((x2[0] - x2[1]) ** 2 + (y2[0] - y2[1]) ** 2)
pixels = np.array([np.average([r1,r2]), np.std([r1,r2])])
print(pixels)
measured distance from Van, from laser center to rounded edge = 4.878 inches above. Alignment stand was 20.85 meters away
angular distance of features is
In [27]:
np.rad2deg(np.arctan(0.1239/20.85))*60*60
Out[27]:
In [28]:
1225/366.78
Out[28]:
confirming with other direction
In [29]:
best_img.imshow()
x1 = [752, 961]
y1 = [756, 756]
x2 = [752, 961]
y2 = [756, 756]
plt.plot(x1, y1)
plt.plot(x2, y2)
Out[29]:
In [30]:
r1 = np.sqrt((x1[0] - x1[1]) ** 2 + (y1[0] - y1[1]) ** 2)
r2 = np.sqrt((x2[0] - x2[1]) ** 2 + (y2[0] - y2[1]) ** 2)
pixels = np.array([np.average([r1,r2]), np.std([r1,r2])])
print(pixels)
perpendicular measured distance from laser center to right edge = 3.15 inches.
angular distance is then
In [31]:
np.rad2deg(np.arctan(0.08/20.85))*60*60
Out[31]:
The plate scale is then
In [32]:
791.4/209
Out[32]:
This seems much too large compared to the other measurement
In [33]:
fits_files
Out[33]:
Back of payload moved UP
In [34]:
up = [dir + 'FOXSI_SAAS_141109_235421.fits', dir + 'FOXSI_SAAS_141109_235428.fits']
saas_up = [image(up[0]), image(up[1])]
Back of the payload moved less UP
In [35]:
less_up = dir + 'FOXSI_SAAS_141110_000232.fits'
saas_less_up = image(less_up)
Back of the payload moved to the RIGHT (looking sunward)
In [36]:
right = dir + 'FOXSI_SAAS_141110_001342.fits'
saas_right = image(right)
In [37]:
cal = dir + 'FOXSI_SAAS_141109_205930.fits'
In [38]:
cal_saas = image(cal)
In [39]:
cal_saas.roi_reset()
plt.imshow(cal_saas.data, origin='upper', cmap=cm.Greys_r)
plt.contour(saas_right.data)
plt.contour(saas_up[0].data)
plt.contour(saas_less_up.data)
Out[39]:
In [40]:
plt.imshow(saas_up[0].data, origin='upper')
Out[40]:
using origin='lower' shows the uncorrected image as shown on the uncorrected SAAS display
In [41]:
fits_files1 = glob.glob('/Users/schriste/Desktop/SAAS/FOXSI_SAAS_141113*')
In [42]:
fits_files1
Out[42]:
In [43]:
sparcs = image(fits_files1[-1])
no_sparcs = image(fits_files1[-4])
In [44]:
sparcs.roi_reset()
plt.figure()
plt.imshow(sparcs.data, cmap=cm.Greys_r, origin='upper')
xy = np.shape(no_sparcs.data)
plt.contour(no_sparcs.data)
Out[44]:
In [45]:
def find_center(img, sigma=0.9, num_circles=1):
edges = filter.canny(img, sigma=sigma)
hough_radii = np.arange(200, 500, 5)
hough_res = hough_circle(edges, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
peaks = peak_local_max(h, num_peaks=2)
if peaks != []:
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius, radius])
best_centers = []
best_radii = []
best_x = []
best_y = []
number_of_best_circles = num_circles
for idx in np.argsort(accums)[::-1][:number_of_best_circles]:
center_x, center_y = centers[idx]
best_x.append(center_x)
best_y.append(center_y)
best_centers.append(centers[idx])
best_radii.append(radii[idx])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 6))
ax1.imshow(edges)
ax2.imshow(img, cmap=cm.gray)
for center, radius in zip(best_centers, best_radii):
circle = plt.Circle((center[1], center[0]), radius, color='r', fill=False)
ax2.add_patch(circle)
plt.axvline(center[1])
plt.axhline(center[0])
print(radius)
print("Calibrated Center X = %s +/- %s" % (np.average(best_x), np.std(best_x)))
print("Calibrated Center Y = %s +/- %s" % (np.average(best_y), np.std(best_y)))
return np.array([[np.average(best_x), np.std(best_x)], [np.average(best_y), np.std(best_y)]])
In [46]:
s1 = find_center(sparcs.data)
In [47]:
s2 = find_center(no_sparcs.data)
sun is 290 pixels across which is equivalent to 30 arcminutes
In [48]:
30*60/(2*290.)
Out[48]:
In [49]:
print(s1)
print(s2)
In [50]:
((s1[0][0] - s2[0][0])*3.1)/60.0
Out[50]:
In [51]:
((s1[1][0] - s2[1][0])*3.1)/60.0
Out[51]:
In [52]:
plt.figure()
plt.imshow(sparcs.data, cmap=cm.Greys_r, origin='upper')
xy = np.shape(no_sparcs.data)
plt.contour(no_sparcs.data)
plt.axhline(s1[0][0], color='r')
plt.axhline(s2[0][0], color='b')
plt.axvline(s1[1][0], color='r')
plt.axvline(s2[1][0], color='b')
Out[52]:
In [53]:
f = dir + 'FOXSI_SAAS_141113_224732.fits'
img = image(f)
In [54]:
img.roi_reset()
fig = plt.figure()
ax = plt.imshow(img.data, cmap=cm.Greys_r, origin='lower')
circle1 = plt.Circle(img.calibrated_center, 15*60/3.1, fill=False, color='r')
circle2 = plt.Circle(img.calibrated_center, 7.5*60/3.1, fill=False, color='r')
circle3 = plt.Circle(img.calibrated_center, 3*60/3.1, fill=False, color='r')
ax.axes.add_artist(circle1)
ax.axes.add_artist(circle2)
ax.axes.add_artist(circle3)
plt.axvline(img.calibrated_center[0])
plt.axhline(img.calibrated_center[1])
Out[54]:
In [55]:
plt.hist(np.ravel(sparcs.data))
Out[55]:
Nov 14-2014
In [56]:
sparcs_on = dir + 'FOXSI_SAAS_141114_174035.fits'
In [57]:
sparcs_on = image(sparcs_on)
In [58]:
#x = find_center(sparcs_on.data)
In [59]:
ax = plt.imshow(sparcs_on.data, origin="upper")
circle = plt.Circle((1296/2, 966/2), 0.5*0.5*60*60/3.1, fill=False)
ax.axes.add_patch(circle)
#plt.plot(x[1][0], x[0][0], 'x')
Out[59]:
In [60]:
np.shape(sparcs_on.data)
Out[60]:
arcminute shift in y direction is
In [61]:
#((x[0][0] - 966/2)*3.1)/60.0
arcminute shift in x direction is
In [62]:
#((x[1][0] - 1296/2)*3.1)/60.0
This agres very well with SPARCS offsets when we centered using the SAAS
In [63]:
filename = dir + 'FOXSI_SAAS_141115_162920.fits'
In [64]:
saas_cal = image(filename)
In [65]:
plt.imshow(saas_cal.data, origin="upper")
Out[65]:
In [66]:
saas_cal.imshow()
saas_cal.overlay()
In [67]:
saas_cal.max_index
Out[67]:
In [68]:
saas_cal.data.max()
Out[68]:
In [69]:
saas_cal.header
Out[69]:
In [70]:
from saaspy.image import image
In [71]:
s = image(filename = dir + 'FOXSI_SAAS_141115_162920.fits')
In [72]:
s.gain_analog
Out[72]:
In [73]:
s.header.get("GAIN_PRE")
Out[73]:
In [74]:
from saaspy.collection import collection
In [75]:
d = collection(fits_files[2:])
In [76]:
temp_files = ['FOXSI_SAAS_141115_235253', 'FOXSI_SAAS_141116_000435', 'FOXSI_SAAS_141116_001714',
'FOXSI_SAAS_141116_002646', 'FOXSI_SAAS_141116_003259', 'FOXSI_SAAS_141116_003901',
'FOXSI_SAAS_141116_004555', 'FOXSI_SAAS_141116_005338', 'FOXSI_SAAS_141116_010201',
'FOXSI_SAAS_141116_011211', 'FOXSI_SAAS_141116_012352', 'FOXSI_SAAS_141116_013642',
'FOXSI_SAAS_141116_014001', 'FOXSI_SAAS_141116_015814']
temp_files = [dir + f + '.fits' for f in temp_files]
temperatures = np.array([-18.9, -35, -35, -30, -25, -20, -15, -10, -5, 0, 5, 9.5, 10.7, 15])
In [77]:
temp_files
Out[77]:
In [78]:
s = []
for f in temp_files:
s.append(image(f))
In [79]:
cl = collection(temp_files)
In [80]:
cl.data.plot(subplots=True)
Out[80]:
In [81]:
x = cl.data['mindex_x']
y = cl.data['mindex_y']
In [82]:
x.min()
Out[82]:
In [83]:
fig, ax = plt.subplots(nrows=3, sharex=True)
ax[0].plot(temperatures, x)
ax[0].set_xlabel('Temperature [C]')
ax[0].set_ylabel('X pixel')
ax[1].plot(temperatures, y)
ax[1].set_ylabel('Y pixel')
ax[2].plot(temperatures, np.sqrt(x ** 2 + y ** 2))
ax[2].set_ylabel('R pixel')
Out[83]:
In [84]:
fig, ax = plt.subplots(nrows=3, sharex=True)
ax[0].plot(temperatures, (x - x.min()) * 3.1, 'o')
ax[0].set_ylabel('X [arcsec]')
ax[1].plot(temperatures, (y - y.min()) * 3.1, 'o')
ax[1].set_ylabel('Y [arcsec]')
ax[2].plot(temperatures, np.sqrt(((x - x.min())* 3.1) ** 2 + ((y - y.min())* 3.1) ** 2), 'o')
ax[2].set_ylabel('R [arcsec]')
fit = np.polyfit( temperatures, np.sqrt(((x - x.min())* 3.1) ** 2 + ((y - y.min())* 3.1) ** 2), deg=1)
yfit = lambda x: fit[1] + x * fit[0]
ax[2].plot(temperatures, yfit(temperatures), 'r-')
ax[2].set_xlabel('Temperature [C]')
print('fit result: y = mx + b, m = %f, b = %f' % (fit[0], fit[1]))
In [85]:
fits_files = glob.glob(dir + 'FOXSI_SAAS*')
In [86]:
print('Number of file = %i' % len(fits_files))
In [87]:
cl = collection(fits_files[50:])
In [88]:
fits_files[0]
Out[88]:
In [89]:
for i, f in enumerate(fits_files):
print(i, f)
image(f)
In [ ]:
cl.data
In [ ]:
cl.data['std'].plot()
In [90]:
files = glob.glob(dir + 'FOXSI_SAAS_141114_1*')
In [91]:
cl = collection(files)
In [92]:
cl.data
Out[92]:
In [93]:
cl.data['max'].max()
Out[93]:
In [94]:
files[0]
Out[94]:
In [95]:
s = image(files[0])
In [96]:
s.imshow()
In [97]:
s.gain_analog, s.gain_preamp, s.exposure
Out[97]:
In [98]:
plt.figure(figsize=(18,4))
ax = plt.hist(s.data.ravel(), bins=np.arange(1,255,1), log=True)
plt.xlim(0,255)
plt.show()
In [105]:
plt.plot(s.data[450, :])
Out[105]:
In [131]:
s = image(files[8])
plt.figure(figsize=(15,15))
plt.imshow(s.data[150:750,300:1000], cmap=plt.cm.Greys_r)
plt.colorbar()
plt.title("Active Regions")
print(s.gain_analog, s.gain_preamp, s.exposure)
print(s.max)
plt.show()
In [101]:
file = dir + 'FOXSI_SAAS_141115_213723.fits'
s = image(file)
In [102]:
s.imshow()
plt.xlim(600, 700)
plt.ylim(400, 500)
plt.plot(645, 453, 'x')
Out[102]:
In [104]:
file = dir + 'FOXSI_SAAS_141115_235253.fits'
s = image(file)
In [105]:
s.imshow()
plt.xlim(600, 700)
plt.ylim(400, 500)
plt.plot(645, 453, 'x')
Out[105]:
In [11]:
s.center_of_mass + np.array([s.roi[2], s.roi[0]])
Out[11]:
In [9]:
s.max_index
Out[9]:
In [10]:
s.roi
Out[10]:
In [4]:
files = glob.glob(dir + "FOXSI_SAAS_141120_09*")
In [5]:
files
Out[5]:
In [6]:
cl = collection(files)
In [7]:
cl.data
Out[7]:
In [8]:
xcom = [x[0] for x in cl.data['com']]
ycom = [x[1] for x in cl.data['com']]
In [9]:
plt.plot(xcom, ycom, 'x')
plt.title("Center of Mass")
Out[9]:
In [10]:
s = image(files[0])
In [11]:
s.calibrated_center
Out[11]:
In [12]:
xcal = s.calibrated_center[1]
ycal = s.calibrated_center[0]
In [13]:
offset_x = [xcal - x for x in xcom]
offset_y = [ycal - x for x in ycom]
In [14]:
plt.plot(offset_x, offset_y, 'x')
plt.plot(np.average(offset_x), np.average(offset_y), 'o')
plt.errorbar(np.average(offset_x), np.average(offset_y), yerr=np.std(offset_y), xerr=np.std(offset_x))
plt.title("Offset in Pixels")
plt.xlabel('X')
plt.ylabel('Y')
Out[14]:
In [15]:
x = []
y = []
for f in files:
xy = find_center(image(f), num_circles=1, sigma=1.3)
x.append(xy[0][0])
y.append(xy[1][0])
print(x)
In [16]:
find_center(image(files[0]))
Out[16]:
In [22]:
plate_scape = 970.38/(290)
print("Plate Scale %f = " % plate_scape)
In [30]:
plt.plot(np.array(offset_x) * plate_scape/60., np.array(offset_y)* plate_scape/60., 'x')
label = "[%f, %f]" % (np.average(offset_x)* plate_scape/60., np.average(offset_y)* plate_scape/60.)
plt.plot(np.average(offset_x)* plate_scape/60., np.average(offset_y)* plate_scape/60., 'o', label=label)
plt.errorbar(np.average(offset_x)* plate_scape/60., np.average(offset_y)* plate_scape/60., yerr=np.std(offset_y)* plate_scape/60., xerr=np.std(offset_x)* plate_scape/60.)
plt.title("Offset in Pixels")
plt.xlabel('X offset [arcmin]')
plt.ylabel('Y offset [arcmin]')
plt.legend()
Out[30]:
In [6]:
file = dir + 'FOXSI_SAAS_141130_110145.fits'
In [7]:
s = image(file)
In [8]:
previous_calibrated_center = np.array([645, 453])
In [9]:
s.max
Out[9]:
In [10]:
s.imshow()
plt.xlim(600,750)
plt.ylim(400, 550)
plt.plot(previous_calibrated_center[0], previous_calibrated_center[1], 'x', label="pre-vibe")
new_calib_center = np.array([676, 472])
plt.plot(new_calib_center[0], new_calib_center[1], 'o', label="post-vibe")
plt.legend()
Out[10]:
In [11]:
print("Pre-vibe calib center", previous_calibrated_center)
print("Post-vibe calib center", new_calib_center)
print("Difference", (previous_calibrated_center - new_calib_center)*s.plate_scale/60.)
print("Shift", np.sqrt(np.sum(np.power(previous_calibrated_center - new_calib_center, 2)))*s.plate_scale/60.0)
In [37]:
file = dir + 'FOXSI_SAAS_141201_065826.fits'
In [38]:
s = image(file)
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], '+', label="post-vibe calib center")
plt.plot(previous_calibrated_center[0], previous_calibrated_center[1], 'x', label="pre-vibe")
plt.legend()
Out[38]:
In [39]:
r = find_center(s, plot=True, num_circles=1)
print(r)
In [40]:
y = r[0][0]
x = r[0][1]
print(x, y)
In [41]:
sun_center = r[0][0:2]
print(sun_center)
In [44]:
s = image(file)
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], '+', label="post-vibe calib center")
plt.plot(x, y, 'x', label='sun center')
plt.plot(previous_calibrated_center[0], previous_calibrated_center[1], '+', label="pre-vibe")
plt.legend()
Out[44]:
shift is compared to pre-vibe alignment (this value given to Chris to make shims)
In [47]:
(np.array([x,y]) - previous_calibrated_center) * s.plate_scale/60.0
Out[47]:
shift is compared to post-vibe alignment (what we really needed to shim)
In [45]:
(np.array([x,y]) - new_calib_center) * s.plate_scale/60.0
Out[45]:
Post-shim image (maybe)
In [ ]:
file = dir + 'FOXSI_SAAS_141201_085533.fits'
s = image(file)
In [50]:
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], '+', label="post-vibe calib center")
plt.plot(previous_calibrated_center[0], previous_calibrated_center[1], '+', label="pre-vibe")
plt.plot(s.center_of_mass[1], s.center_of_mass[0], 'x', label='sun center')
plt.legend()
Out[50]:
In [27]:
file = dir + 'FOXSI_SAAS_141202_072511.fits'
s = image(file)
In [28]:
r = find_center(s, plot=True, num_circles=1)
sun_center = r[0][0:2]
In [30]:
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], '+', label="post-vibe calib center")
plt.plot(s.center_of_mass[1], s.center_of_mass[0], '+', label='center of mass')
plt.plot(sun_center[1], sun_center[0], 'xb', label='sun center')
plt.legend()
Out[30]:
In [51]:
file = dir + 'FOXSI_SAAS_141202_074421.fits'
s = image(file)
In [52]:
fig, ax1 = plt.subplots(1, 1, figsize=(8, 6))
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], 'o', label="post-vibe calib center")
#plt.plot(s.center_of_mass[1], s.center_of_mass[0], 'x')
#plt.plot(sun_center[1], sun_center[0], 'xb', label='sun center')
sun_center = (620, 485)
circle = plt.Circle(sun_center, radius=280, color='r', fill=False)
plt.plot(sun_center[0], sun_center[1], 'x', label='Sun Center')
ax1.add_patch(circle)
plt.legend()
plt.show()
In [53]:
(new_calib_center - sun_center) * s.plate_scale/60.0
Out[53]:
checking next image which should be the same
In [54]:
file = dir + 'FOXSI_SAAS_141202_074414.fits'
s = image(file)
In [55]:
fig, ax1 = plt.subplots(1, 1, figsize=(8, 6))
s.imshow()
plt.plot(new_calib_center[0], new_calib_center[1], 'o', label="post-vibe calib center")
#plt.plot(s.center_of_mass[1], s.center_of_mass[0], 'x')
#plt.plot(sun_center[1], sun_center[0], 'xb', label='sun center')
sun_center = (620, 485)
circle = plt.Circle(sun_center, radius=280, color='r', fill=False)
plt.plot(sun_center[0], sun_center[1], 'x', label='Sun Center')
ax1.add_patch(circle)
plt.legend()
plt.show()
In [ ]: